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COORDINATES  DERIVED  FROM  SPECIFIED  GEOGRAPHIC  COORDINATES 
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Introduction:  For  many  purposes  it  is  useful  to  have  a  fast  means  of  estim¬ 
ating  directional  cosmic  ray  cutoffs  pertaining  to  any  specified  direction 
and  location.  The  derivation  of  "real**  geomagnetic  field  cutoffs  by  means 
of  computer  calculations  (see  McCracken  et  al.  1962;  Freon  and  McCracken, 
1962;  Shea  et  al,  1965;  Cooke  and  Humble,  1979;  for  example)  is  a  slow  op¬ 
eration,  and  expensive  in  computer  time.  The  Stormer  cutoff  function  (Stormer 
(1930,  1955)),  which  expresses  the  dependence  of  the  Stormer  cutoff  on  location 
and  direction  in  a  dipole  approximation  to  the  earth's  field,  offers  a  means 
of  determining  cutoffs  which  is  many  orders  of  magnitude  faster  than  the 
trajectory  tracing  method,  and  one  which  is  sufficiently  precise  to  produce 
useful  estimates  of  cutoffs  in  non-critlcal  applications. 

The  fundamental  imprecision  in  the  Stormer  expression  in  representing 
the  real  field  cutoffs,  in  particular  due  to  failure  to  take  into  account 
higher  order  field  harmonic  terms,  or  the  width  and  transparency  of  the 
penumbra,  is  normally  exacerbated  by  the  use  of  earth  centered  geographic  or 
geomagnetic  coordinates  when  invoking  the  expression  (because  the  earth's 
equivalent  dipole  is  not  earth  centered).  It  is  possible  to  appreciably 
improve  the  accuracy  of  the  cutoff  estimates  by  using  "magnetic”  coordinates 
(l.e.  offset  dipole  latitude,  longitude,  zenith,  and  azimuth)  when  employing 
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the  expression,  and  in  this  way  to  take  into  account  the  offset  and  tilt  of 
the  earth's  equivalent  dipole.  By  this  means,  inherently,  the  effect  of 
Ignoring  the  higher  order  field  terms  is  minimized.  Smart  and  Shea  (1977) 
have  discussed  the  advantages  of  using  offset  dipole  coordinates  In  conjunc¬ 
tion  with  the  Stormer  expression,  and  the  use  of  this  expression  for  inter¬ 
polating  cosmic  ray  cutoffs  over  Intervals  within  which  precise  calculated 
values  do  not  exist. 

A  self  contained  system  for  transforming  from  geographic  coordinates  to 
offset  dipole  coordinates  is  described  in  this  report.  A  printed  copy  of  the 
computer  program  containing  this  procedure,  which  uses  the  calculated  coor¬ 
dinates  to  derive  cutoff  estimates  from  the  Stormer  expression,  is  appended 
to  the  report. 

Discussion;  The  coordinate  transformation  has  a  number  of  stages,  which  are 
Individually  described  in  the  following: 

1)  A  coordinate  conversion  Is  used  to  determine  the  offset  dipole  lati¬ 
tude,  longitude,  and  radius  from  the  nominated  geographic  latitude  and  longi¬ 
tude,  and  geocentric  altitude.  This  conversion  takes  into  account  the  offset 
and  inclination  of  the  earth's  equivalent  dipole  for  any  required  epoch,  and 
assumes  that  the  earth  is  an  oblate  spheroid  of  eccentricity  0.00674.  The 
angle  conversion  equations  are  as  follows: 

offset  dipole  longitude  $  ”  arctan  *1°  jt  COB  ~  y  co8  j>) 

C 

offset  dipole  latitude  9  -  arctan  (cos  $  tan  a)  I 

C 

offset  dipole  radius  R'  ” 

cos  0  cos  $ 


where  a  *  arctan  (F/G)  +  a 
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C  a  G  C08a 
cos(  a  -  a) 

and  F  »  R  sin  X  -  x 

G  ■  y  sinb  +  R  cos  X  cos  ( —  c) 

X  and  t|>  are  the  geographic  latitude  and  longitude  respectively;  R  Is  the 
geographic  radius  at  the  specified  location;  x  Is  the  displacement  of  the 
equivalent  earth  dipole  above  the  equatorial  plane;  y  is  the  displacement  of 
the  dipole  from  the  center  of  the  earth  In  a  direction  parallel  with  the 
geographic  equatorial  plane;  a  is  the  Inclination  of  the  dipole  axis  from  a 
direction  parallel  to  the  geographic  N-S  axis;  b  is  related  to  the  angle 
between  the  zero  magnetic  longitude  and  the  geographic  longitudinal  direction 
towards  which  the  dipole  is  displaced;  and  c  is  the  geographic  longitudinal 
direction  parallel  to  which  the  direction  zero  offset  dipole  longitude  lies. 

A  procedure  for  establishing  the  position  of  the  equivalent  dipole  with 
respect  to  the  geocentric  coordinate  system  at  any  epoch  is  presented  by  Smart 
and  Shea  (1977).  This  procedure,  which  makes  use  of  the  low  order  terms  in 
the  spherical  harmonic  representation  of  the  geomagnetic  field,  is  discussed 
more  fully  by  Roederer  (1972),  and  Akasofu  and  Chapman  (1975).  The  values  of 
x,  y,  a,  b  and  c  used  in  the  presently  reported  coordinate  transformation  can 
thus  be  determined  as  required.  The  equivalent  dipole  position  determination 
is  included  as  a  subroutine  in  the  computer  program  appended  to  this  report. 

2)  The  following  steps  rely  on  the  use  of  a  particular  means  of  specifying 
the  relative  position  of  two  points  on  the  surface  of  a  sphere.  In  particular, 
two  angles  are  used,  one  (o)  defines  the  angle  between  the  great  circle  con¬ 
necting  the  two  points  and  the  meridian  line  intersecting  one  of  them,  and 
the  other  (y)  is  the  angle  subtended  at  the  center  of  the  sphere  by  the  two 
points  (see  figure  1). 


* 
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Figure  l.  Diagram  defining  the  angles 
a  and  Y  used  to  express  the 
relative  position  of  the  two 
points  Pi  and  P£  on  a  spheri¬ 
cal  surface.  C  is  the  centre 
of  the  sphere. 


If  the  two  points  are,  as  specified  by  latitude  and  longitude,  (Xi,i|>i) 
and  cllen  "  Y  are  given  by 


a  m  arctan  (_ 


sin  ( \1>  -  \1»  )  cos  X 
2  1  2 


sin  X  cos  X  -  sin  X  cos  X  cos  ( -  \li  ) 
12  12  2  1 


Y  ■  arcos  (cos  (f2  ~  tl)  cos  Xi  cos  X2  +  sin  Xi  sin  X2) 


J 


II 


Conversely,  if  the  location  of  one  point  is  known  (Xi,\|»i,  say),  then  it  is 
possible  to  determine  the  latitude  and  longitude  of  a  second  whose  position 
is  defined  in  terms  of  a  specified  o  and  y,  by  means  of 


latitude  ■  arcsin  (cos  a  sin  y  cos  Xj  +  cos  y  sin  Xi) 


cos  y  "  sin  X  sin  latitude 
longitude  *  1  +  arcos  ( _ 1 _ ) 

cos  X  cos  latitude 
1 


III 


3)  The  calculation  of  the  magnetic  zenith  and  azimuth  uses  the  strategy 
of  setting  up  a  vector  of  known  length  (R")  pointing  in  the  direction  of 
interest.  The  position  the  tall  of  the  vector,  being  the  location  of  inter¬ 
est,  is  of  course  known  in  both  geographic  and  magnetic  coordinates  (the 
latter  determined  by  step  1),  and  similarly  both  geographic  and  offset  dipole 
radius  values  pertaining  to  this  point  are  known. 
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The  geographic  coordinates  of  the  position  of  the  head  of  the  vector  can 
be  determined  by  using  calculated  o  and  y  values,  derived  as  follows  (which 
values  pertain  to  the  projection  of  the  vector  onto  a  spherical  surface). 


y  m  arctan  (R"  sin  ze/(R  +  R”  cos  ze)) 


IV 


o  -  az 

Having  evaluated  these  angles,  the  geographic  latitude  and  longitude  of  the 
vector  head  can  be  calculated  using  the  relationships  III.  The  offset  dipole 
latitude  and  longitude,  and  the  distance  R*  of  the  vector  head  from  the  dipole 
are  then  calculated  using  relationships  I.  Now,  having  the  offset  dipole 
coordinates  of  both  ends  of  the  vector,  the  angles  a  and  y  relating  the  two 
points  (in  the  offset  dipole  frame  of  reference)  can  be  calculated  by  means 
of  the  relationships  II.  The  azimuth  of  the  vector  in  the  offset  dipole 
frame  of  reference  is  simply  the  a  value,  whilst  the  zenith  angle  is  given  by 


ze  -  arcos  (R*  cos  y  -  R* ) 


where  R*  is  the  offset  dipole  radius  of  the  vector  tail  (i.e.  the  distance  of 
the  vector  tail  from  the  offset  dipole  center). 


Conclusion:  By  using  the  coordinate  transformation  procedure  described,  the 
position  (latitude,  longitude,  and  radius)  of  the  site  location  relative  to 
the  offset  dipole,  and  the  zenith  and  azimuth  pertaining  to  the  direction  of 
Interest,  may  be  calculated  from  the  specified  geographic  coordinates  and 
geocentric  altitude.  At  each  step  during  the  derivation  of  these  parameters 
tests  are  performed  to  ensure  that  calculated  angle  values  lie  in  the  correct 
quadrant,  and  appropriate  corrections  are  made  if  they  do  not.  Having  thus 
determined  the  angles  and  distance  relative  to  the  equivalent  dipole  centered 
frame  of  reference  the  Stormer  equation  can  be  invoked  with  greatest  possible 
precision. 


The  appropriately  normalized  Stormer  expression  is  as  follows: 
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cutoff  ”  59.4  cos  8 

R' ^  (l  +  /  1  -  cob^  9  sin  az  sin  ze)2 

This  expression  takes  in  offset  dipole  latitude  (9),  zenith  (ze),  azimuth  (az, 
measured  clockwise  from  magnetic  north),  and  radial  distance  from  the  effec¬ 
tive  dipole  center  (R’);  and  produces  cutoff  rigidity  values,  in  units  of  GV. 

The  appendix  contains  a  printed  copy  of  the  computer  program  which  carries 
out,  in  the  way  described,  the  entire  task  of  producing  a  Stormer  cutoff  esti¬ 
mate  from  the  specified  input  geographic  position  and  direction  parameters. 
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PROGRAM  STORM 

DOUBLE  PRECISION  ZE,  AZ,  LAT, LON. ZER. AZR.  LATR. LONR.  ZERO.  DTOR 
DOUBLE  PRECISION  MLATA.  MLONA,  ZEMA. AZMA, ZENEW. AZNEW. GAM. SIG 
DOUBLE  PRECISION  DSIN,  DCOS. DTAN. DASIN. DACOS. DATAN.  DABS,  DSGRT 
DOUBLE  PRECISION  DFLOAT.  INCR.  ERAD, DIFF, A, B. X, Y, R, RDEL. ALT 
DOUBLE  PRECISION  GAMVAL.  MAOLAT, MACLON. RMAG, CUTOFF. PI,  P ID2.  PIB2 
COMMON  /PIANO/PI,  PID2, PIB2, ZERO 
COMMON  /ANCLES/  SIG, GAM 
COMMON  /POS/  X, Y, A, B, DIFF, ERAD 

COMMON  /MACIO/  ZENEW,  AZNEW.  R,  RMAG. CUTOFF. MAGLAT. MACLON, GAMVAL. 

X  RDEL. ALT 

C  ANGLE  CONVERSION  FACTORS  AND  CONSTANTS  ARE  PRESET  IN  THE  FOLLOWING 
PI  *  3  141 59265359D *-00 
PID2  =  PI/2.  OD+OO 
P1B2  -  P 1*2.  OD+OO 
DTOR  »  PI/ 180.  OD+OO 
RDEL  *  100  OD+OO 
I NCR  *  1  0D-10 
ZERO  *  0  OD+OO 

C  POSITION  AND  INCLINATION  OF  OFFSET  DIPOLE  NOW  DEDUCED 
CALL  POSITN 

C  IF  DESIRED,  INSERT  OTHER  VALUES  OF  A,  B,  X,  Y,  DIFF  AND  ERAD  .Hi. 

A  *  A*DTOR 
B  =  ii+DTOR 
DIFF  *  DiFF*DTOR 
C  INPUT  DATA  ACCEPTED  HERE. 

100  TV PE  1000 

ACCEPT  *. LAT, LON, ZE. AZ. ALT 

C  TERMINATE  RUN  IF  LATITUDE  VALUE  ENTERED  IS  OUTSIDE  THE  POSSIBLE  RANGE 
IF  (DABS (LAT >  .  GT.  90.  OD+OO)  CO  TO  200 
ZER  =  ZE*DTOR+INCR 
AZR  =  AZ*DTOR+INCR 
CALL  RANGE (AZR) 

LATR  *  i_AT*DTOR+INC.R 
LONR  =  LON+DTOR+INCR 

C  NOW  DEDUCE  OFFSET  DIPOLE  COORDINATES  OF  DIRECTION. 

CALL  MAGCO(LATR,  LONR, ZER,  AZR) 

C  CALCULATED  ANGLES  ARE  NOW  CONVERTED  FROM  RADIANS  TO  DEGREES 
MLATA  =  MACLAT/DTOR 
MLONA  *  MAGLON/DTOR 
ZEMA  =  ZENEW/DTOR 
AZMA  *  AZNEW/DTOR 

C  RESULTS  ARE  PRINTED  OUT  BY  THE  FOLLOWING  INSTRUCTIONS. 

WRITE <6. 2000) 

WRITE! 6, 4000) 

WRITE (6.  5000)  LAT,  LON,  ZE,  AZ,  ALT 

WRITE <6, 6000)  MLATA  MLONA, ZEMA. AZMA. RMAG, CUTOFF 

WRITE  < 6, 3000) 

GO  TO  100 

1000  FORMAT ('  ENTER  LAT,  LON,  ZE,  AZ,  AND  ALT.  ',*) 

2000  FORMAT (IX) 

3000  FORMAT (IX,  /) 

4000  FORMAT  ( 17X, 45H  LATITUDE  LONGITUDE  ZENITH  AZIMUTH. 

X  36H  ALTITUDE  DIP.  RAD.  CUTOFF) 

5000  FORMAT  ( 13H  GEOGRAPHIC:  , 5F12.  2) 

6000  FORMAT  ' 13H  OFFSET  DIP  :  . 4F12.  2,  12X, 2F12.  2) 

200  STOP 
END 
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.'0710 
00720 
00730 
00740 
00750 
00760 
00770 
00780 
00790 
■'0800 
'0810 
00820 
00830 
00840 
00850 
00860 
00870 
00880 
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COP>00 
r’0910 
00920 
00^30 
009  40 
00950 
00960 
00970 
00980 
00990 

:  iooo 
01010 
01020 
01030 
01040 
01050 
0 1 060 
C  1 070 
C  1 080 
01090 
O'  1  1 00 
OHIO 
1  1 20 
0 1 1 30 
:  i  i  40 
0  1  1 50 
01160 
01170 
0  1  1  90 
01190 
C  1200 
01210 
01220 
01230 
11240 
01250 
01260 
01270 
01280 
01290 
01300 
01310 


SUBROUTINE  RANGE  < PARAM) 

C 

C  THIS  SUBROUTINE  BRINGS  ANGLE  VALUES  TO  WITHIN  0-360  DECREE  RANGE 
C 

COMMON  /PIANO/PI. PID2,  PIB2, ZERO 
DOUBLE  PRECISION  PARAM,  PI , PID2, R1B2, ZERO 
IF  < PARAM  .  GE.  PIB2)  PARAM  ■  PARAM-PIB2 
IP  (PARAM  .  LT.  ZERO)  PARAM  -  PARAM+PIB2 

RETURN 

END 


SUBROUTINE  CHECK  < PARAM ) 

C 

C  THIS  SUBROUTINE  PREVENTS  EFFECT  OF  ROUND-OFF  ERRORS  ETC  CAUSING 
C  THE  ARGUMENT  OF  AN  ANGLE  FUNCTION  TO  EXCEED  1  00  MAGNITUDE 
C 

DOUBLE  PRECISION  PARAM 

IF  < PaRaM  .  GT.  +1  OD+OO)  PARAM  *  1  CD+OO 
IF  (PARAM  LT  -1  OD+OO)  PARAM  *  -l.OB+CO 
RETURN 
END 


SUBROUTINE  CONV ( LAT 1 , LAT2, L0N1 . L0N2 ) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  GAMMA  AND  SIGMA  ANGLES  RELATING  "rHE 
C  Two  Earth-  OR  DIPOLE -CENTERED  VECTORS  DIRECTED  TOWARDS  (LAT1, L0N1 ) 

C  AND  '■ '  .AT2  ■  L0N2  ) . 

C 

DCUBlE  PRECISION  ChvAL.  GAM, SIG. DS1N, DASIN- DCOS, DACOO,  PI, PID2, PIB2 
DOUBLE  PRECISION  LAT 1 ,  LAT2,  L0N1 >  L0N2- VAL, BTLR,  GAMD, SIGD,  VALD,  ZERO 
C OMMON  /P I ANG/  PI, PID2. PI B2, ZERO 
COMMON  /ANGLES/  SIG, GAM 

CHVAL  *  DCOS <L0N1-L0N2)*DC0S( LAT 1 ) *DCGS ( LAT2) +DSIN ( LAT 1 /* 

X  DSIN1LAT2) 

CALL  CHECK (CHVAL I 
GAM  -  DACOS< CHVAL' 

SIG  =  ZERO 

IF  (Gam  EQ  ZERO;  CO  TO  100 

CHVAL  =  D5IN(L0N2-LGN1 )*DC0S(LAT2> /DSIN(GAM) 

CALL  CHECK (CHVAL) 

SIG  -  DASIN <, CHVAL ) 

100  CONTINUE 

CHVAL  *  DCOS  (  GAM  .)  *-D5IN  <  LAT  1  ) 

CALL  CHECK (CHVAL) 

VAL  *  DACOS  <  CHVAL ' 

IF  ((PI D2-LAT2 )  GE  VAL)  SIG  *  PI-SIG 
CALL  RANGE (SIC) 

RETURN 

END 


9 


01320 
01330 
01340 
0 1 350 
01360 
01370 
01380 
01390 
0 1 400 
01410 
01420 
01430 
01440 
01450 
01460 
01470 
01480 
0 1 490 
0 1 500 
01510 
0 1  520 
0 1  530 
0 1  540 
0 1 550 
0 1  560 
01570 
Cl  580 
01590 
0 1  600 
01610 
01620 
01630 
01640 
01650 
01660 
01670 
0 1 680 
01690 
01700 
01710 
01720 
0 1 730 
01740 
01750 
0 1 760 
01770 
01780 
01790 
01800 
01810 
01820 
01830 
01840 
01850 
01860 
01870 
01880 
01890 
01900 
01910 
01920 


SUBROUTINE  CALC <L AT,  LON, CAM, SIC, LAT  ,  , LONN) 

C 

C  THIS  SUBROUTINE  DEDUCES  THE  LATITUDE  AND  LONCITUDE  OF  A  DIRECTION 
C  WHICH  HAS  A  NOMINATED  SIGMA  AND  GAMMA  DEVIATION  FROM  A  SPECIFIED 
C  LATITUDE  AND  LONCITUDE 
C 

DOUBLE  PRECISION  SIS,  CAM,  LAT, LON, LATT, LONN, FUNC.  PI,  PID2,  PIB2.  ZERO 
DOUBLE  PRECISION  DSIN,  DCOS,  DASIN,  DACOS,  DABS,  DFLOAT,  VAL.  CHVAi. 
COMMON  /PIANG/  PI,  PID2,  PIB2, ZERO 

CHVAL  =  DCOS  ( SIC )  *D5IN( GAM )  *DCOS  ( LA  t  )  ♦DCOS (GAM )  *DSIN(LAT ) 

CALL  CHECK! CHVAL) 

LATT  *  DAS I N  <  CHVAL  > 

FUNC  =  1  OD-t-OO 

Ir  (SIG  GT  PI)  FLi-JC  ■  -1  OD+OO 

CHVAL  =  (  DC  OS  i  CAM  >  -DS I N  (  LAT )  *DS  i  N  •.  L A  TT  )  »  /DCOS.  (  LAT  ) /DCOS  (l  ATI  ) 

CALL  CHECK (CHVAL)  \ 

Li.'NN  =  LON+FUNC*DAC  OS  (CHVAL) 

IF  ( LATT  LT  PID2>  CO  TO  100 
LATT  =  PI -LATT 
Ll'NN  =  LONN-PI 

ICO  Ir  (LATT  CT.  -PID2'  CO  TO  200 
LATT  =  -LATT— P I 
LONN  =  LONN-PI 
200  CALL  *ANCE(LONN) 

Rt  TURN 
END 


SUBROUTINE  GETOMA<  LAT,  LON. RVEC  > 

C 

C  THIS  SUBRC'UT  INE  TAKES  IN  GEOGRAPHIC  LATITUDE  fc  LONGITUDE  AND 
C  OUTPUTS  OFFSET  DIPOLE  LATITUDE  8.  LONCITUDE.  LOCAL  GEOGRAPHIC 
C  RADIUS.  AND  DISTANCE  FROM  THE  OFFSET  DIPOLE. 

C 

DOUBLE  PRECISION  P  I  PID2,  P  IB2,  A  B  ■  X  V-  DIFF,  LAT,  LON.  /ENEW,  AZnE**,  ALT 

DOUBLE  PRECISION  R .  RMAG,  CUTOFF. MAGLAT, MAGLON.  GAMVAL,  RVEC.  RDEL-  VAL 

DOUBLE  PRECISION  DSIN, DCOS, DTAN, DASIN- DACOS. DATAN, DSQRT,  DABS.  ZERO 

DOUBLE  PRECISION  ePAD 

COMMON  /PIANG/  PI, PID2, PIB2, ZERO 

COMMON  /POS/  X. Y. A, B, DIFF, ERAD 

COMMON  MAG  10/  Z ENEW ,  AZNEW,  R,  RMAG-  CUTOFF,  MAGLAT,  MAGI  ON,  GAMVAL. 

X  PDEL-  ALT 

R  a  ALT+-6356  7747D+00/DSQRT!  1.  OD-t-CO-O  00673966D+00* 

X  DSIN(PID2-LAT)**2) 

II-  (RVEC  GT  ZERO)  R  *  RVEC 
AVAL  *  R*DSIN(LAT )— X 

VAL  =  Y*DSIN(B)+R*CCOS<LAT>*DCOS(LDN-DIFF> 

GAMVAL  =*  DATAN (AVAL/ VAL )+A 

BVAL  *  VAL*DCOS< GAMVAL) /DCOS (GAMVAL-A) 

MAGLON  «  DATAN ( ( R*D3IN( LON-DIFF > *DCOS < LAT ) -Y*DCOS(B ) )/BVAL) 

IF  (BVAL  LT  ZERO)  MAGLON  *  MACLON+PI 
CALL  RANGE (MAGLON) 

MAGLAT  *  DAT AN (DC OS ( MAGLON >*DT AN < GAMVAL ) ) 

RMAG  *  DABS (BVAL /DCOS (MAGLAT) /DCOS ( MAGLON ) ) 

RETURN 

END 


S''BRO.  TINE  MAGC 0  <  L  AT ,  LON .  ZE. AZ) 


0  1 930 
•1940 
01950 
01 960 
11970 

•  1 980 
>1990 

02000 
02010 
02020 
02030 
02040 
02050 
02060 
02070 
12080 
02090 
02100 
02110 
02 1 20 
■•*>2 1 30 
02140 
02 1 50 
02 1 60 
•02170 
02180 
02190 
02200 
02210 
>12220 
•.•2230 
>-•2240 
12250 
•.•2260 
12270 
12280 
•02290 
>02300 
02310 
02320 
02330 
02340 
02350 
02360 
02370 

•  -•2380 
1>2390 
02400 
02410 
02420 
12430 
02440 
12450 
02460 
02470 
02480 
02490 
12500 
02510 
02520 
02530 


C 

C  THIS  SUBRlUTINE  takes  a  qiven  set  of  geographic  coordinates 
C  <LA>  •  LON.  ZE- AZ  >  AND  DEDUCES  THE  CORRESPONDING  SET  OF  GF^fcT 
C  DIPC'LE  COORDINATES 
C 

DOUBLE  PRECISION  A.  EU  X,  V,  ZENEW  AZNEW, R. RMAG,  CUTOFF 

DOUBLE  PRECISION  ZEM.  AZM.  MAGLAT. MAGI  ON,  GAMVAL,  RDEL,  RV-  CHVAL 

DDUBLE  PRECISION  DS1N, DCOS, DTAN- DAS  IN. DACGS,  DATAN. LAT.  LON-  ZE  aZ 

DOUBLE  PRECISION  MLAT, MLON, EP.  RVEC. GAM. SIG, LATVEC. « QNVEC,  Ai  T 

COMMON  /POS/  X,  Y  A, B, DIFF, ERAD 

COMMON  MACI07  ZENEW. AZNEW, R, RMAG, CUTOFF. 

X  MAGLAT , MAGLON- GAMVAL,  RDEL . ALT 

COMMON  /ANGLES/  SIG. GAM 

C  DEDUCE  OFFSET  DIPOLE  LATITUDE  AND  LONGITUDE  FROM  GEOGRAPHIC  VALUES 
C  ALL  GET OMA < LAT ,  LON  -l.OD+OO) 

<=>V  *  RMAG 
M!  AT  -  MAGLAT 
MUON  =  MAGLON 

C  DEDUCF  GAMMA  AND  SIGMA  APPLYING  TO  SAMPLE  VECTOR 
EPS  =  DaTAN>:RDEL«DSIN(ZE)/<R+RDEL»DCOS<ZEO  ) 

PVEC  =  >:R+RDEL*DCOS<ZE))/DCOS(EPSl 
GAM  ^  EPS 
SIG  =  AZ 

C  CALCULATE  GEOGRAPHIC  LAT  S<  LON  OF  VECTOR  HEAD 
CALL  ' ALC^LAT, LON, GAM,  SIG,  LATVEC. LONVEC ) 

C  CONVERT  THESE  GEOGRAPHIC  COORDINATES  INTO  OFFSET  DIPOLE  COORDINATES 
C  ALL  GET  OMA <  LATVEC  -  LONVEC ,  RVEC ) 

C  CALCUI  ATE  GAM  *  SIG  APPLYING  TO  OFFSET  DIPOLE  COORDINATE  SYSTEM 
CALL  IGNV ( ML AT ,  MAGLAT, ML ON, MAGLON) 

C  SO,  DEDUCE  OFFSET  DIPOLE  ZENITH  AND  AZIMUTH  . 

AZNEW  *  SIG 

CHVAl.  =  ( RMAG*DCOSi GAM )-RV> /RDEL 
CALL  CHEC* (CHVAL) 

CENEW  =  DACQS( CHVAL) 

C  NOW  COMPUTE  CUTOFF  VALUE  FROM  STORMER  EXPRESSION. 

CALL  CL’TCAL ( MLAT ,  MLON.  ZENEW,  AZNEW,  PV,  CUTOFF) 

MAGLAT  =  MLAT 
MAGLON  =  MLON 
FMAG  =  RV 
RETURN 
END 


SUBROUTINE  CUTCALfLAT, LON, ZE, AZ. R. CUTOFF) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  STORMER  CUTOFF  CORRESPONDING  TO  THE 
C  SPECIFIED  OFFSET  DIPOLE  LATITUDE.  LONGITUDE, ZENITH  AND  AZIMUTH 
C 

DOUBLE  PRECISION  LAT, LON, ZE, AZ, P. CUTOFF 

DOUBLE  PREC ISION  DSIN,  DCOS, DTAN, DA3IN, DACOS, DATAN, DSQRT 
CUTOFF  =  59  4D*00*( DCOS < LAT) >**4/< (fi/6400  0D+00>**2* 

X  ll  OD+QO+DSQRT  ( 1 .  0D*00-  ( DCOS  ( LAT  >  )  *#3*DS  I N  ( AZ  >  * 

X  DSIN(ZE) ) >**2> 

RETURN 

END 
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02540 

<-'2550 

32560 

C2570 

02580 

02590 

02600 

02610 

02620 

02630 

02640 

02650 

02660 

02670 

02680 

02690 

02700 

02710 

02720 

02730 

02740 

02750 

02760 

02770 

02780 

32790 

02800 

02810 

02820 

02830 

02840 

02850 

O2860 

-2870 

O2880 

02890 

02900 

02910 

02920 

(->2930 

02940 

02950 

02960 

02970 

02980 

02990 

03000 

03010 

03020 

03030 

03040 

03050 

03060 

33070 

03080 

03090 

03100 

03110 

03120 

03130 


SUBROUTINE  POSITN 
C 

C  THIS  SUBROUTINE  DETERMINES  THE  POSITION  OP  THE  EQUIVALENT  DIPOLE 
C  IN  GEOGRAPHIC  COORDINATES.  AND  RETURNS  THE  DEFINING  PARAMETERS 
C  X,  V,  A,  B.  &  DIFF. 

C 

DOUBLE  PRECISION  X.  V,  A,  B.  DIFF, ERAD. P I . P ID2. P IB2. ZERO 
DOUBLE  PRECISION  GOl.  Qll,  002.  Gl2.  G22,  KOI.  HI  1,  H02.  H12.  H22 
DOUBLE  PRECISION  Cl  1,  PHI,  PHID,  H,  LO.  LI  L2.  E,  DENOM.  XED.  YED.  ZED 
DCU8LE  PREC I SION  DSGRT ,  DATAN.  DFLOAT, CUBERT,  THETA.  THfcTAD. PSI  PSID 
COMMON  POS/  X,  Y- A.  B,  DIFF,  ERAD 
COMMON  /PIANO/  PI.  PID2. PIB2, ZERO 

C  FIE.  D  COEFFICIENTS  FOR  THE  EPOCH  OF  INTEREST  ARE  INSERTED  HERE. 

C  THE  FOLLOW  I NO  COEFFICIENTS  ARE  FOR  THE  IGRF65  FIELD 
EPAD  -  &371.2D+00 
GO 1  *  -0  30339D+00 
Gtl  -  -0.  02123D+00 
S02  -  -0.  01654D+00 
Gl2  *  0  02994D+00 
G22  +  0  01 567D+00 
M01  =  0  OD+OO 
Mil  -  0  05758D+00 
H02  =  0  OD+OO 
HI 2  -  -0.  02006D+00 
H22  =«  0  00130D+00 
Oil  =»  DSQRT(GU**2+H1 1**2) 

FHI  *  DATAN(H11/Gl 1 > 

THETA  =*  DATAN<C  11/001  ) 

H  *  "SORT (001**2+01 1**2+Hl 1**2  > 

CUBEP T  *  3.  OD+OO** ( t .  OD+OO/DFLOAT ( 3 > ) 

lC  *  2  OD+00*G01+GG2+CUBERT* (G1 1+G12+H1  1*H12> 

LI  »  -Gtl*002+CUBERT*<G01*G12+C11*G22+H11*H22> 
l2  »  -HI 1*HC2+CUEERT*(G01*H12-Hl 1*022+01 1*H22> 

E  =  <  L0+G01+L1*G 1 l  +  L2*Hl 1 ) / (4  0D+00+H*H ) 

DENOM  =  3  0D+00*H*H 

*FD  --  (Ll-Gll*E)*ERAD/DENOM 

•FD  -  ( L2-H1 1*E ) +ERAD/DENOM 

ZtD  =  (l0-G01*E)*ERAD/DEN0M 

DTQR  *  PI /ISO  OD+OO 

PSI  =  DATAN (YED/ XED) 

IF  (THETA  OT.  ZERO)  00  TO  100 
THETA  *  -THETA 
PSI  =  PSI+PI 
CALL  RANGE (PS I > 
too  A  «  THETA/DTOR 
DTFF  -  PHI/DTOR 
PSID  -  PSI/DTOR 
B  *  PSI-PHI-PID2 
CALL  RANGE(B) 

CALL  RANGE (PHI) 

E  -  3 /DTQR 
DIFF  *  PHI/DTOR 
X  -  ZED 

Y  -  DSGRT <XED*XED+Y£D*YED> 

C  CONSTANTS  DEFINING  THE  POSITION  AND  OR IFNTAT ION  OF  THE  PQUIVALEUT 
C  DIPOLE  ARE  PRINTED  HERE 

WRITE(6.  1000)  X.Y,  A.B.DIFF 
1000  FORMAT ( 9F12  2) 

RETURN 

END 


